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Abstract — Finding control laws (pulse sequences) that can 
compensate for dispersions in parameters which govern the 
evolution of a quantum system is an important problem in the 
fields of coherent spectroscopy, imaging, and quantum informa- 
tion processing. The use of composite pulse techniques for such 
tasks has a long and widely known history. In this paper, we 
introduce the method of Fourier synthesis control law design 
for compensating dispersions in quantum system dynamics. We 
focus on system models arising in NMR spectroscopy and NMR 
imaging applications. 

I. Introduction 

Many applications in the control of quantum systems 
involve controlling a large ensemble using the same control 
signal [1], [2]. In many practical cases, the elements of the 
ensemble show dispersions or variations in the parameters 
which govern the dynamics of each individual system. For 
example, in magnetic resonance experiments, the spins in 
an ensemble may have large dispersions in their resonance 
frequencies (Larmor dispersion) or in the strength of the 
applied radio frequency fields (rf inhomogeneity) seen by 
each member of the ensemble. Another example is in the 
field of NMR imaging, where a dispersion is intentionally 
introduced in the form of a linear gradient [2], and then 
exploited to successfully image the material under study. 

A canonical problem in the control of quantum ensembles 
is the design of rf fields (control laws) which can simulta- 
neously steer a continuum of systems, characterized by the 
variation in the internal parameters governing the systems, 
from a given initial distribution to a desired final distribution. 
Such control laws are called compensating pulse sequences 
in the Nuclear Magnetic Resonance (NMR) literature. From 
the standpoint of mathematical control theory, the challenge 
is to simultaneously steer a continuum of systems between 
points of interest using the same control signal. Typical 
designs include excitation and inversion pulses in NMR 
spectroscopy and slice selective pulses in NMR imaging [2], 
[5], [6], [7], [8], [9], [10], [11], [12], [13], [14]. In many 
cases, one desires to find a control law that prepares the final 
state as some desired function of the parameters. A premier 
example is the design of slice selective pulse sequences in 
magnetic resonance imaging applications, where spins are 
excited or inverted depending upon their physical position 
in the sample under study [2], [3], [15], [16], [17], [18]. In 
fact, the design of such pulses is a fundamental requisite 
for almost all magnetic resonance imaging techniques. This 
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paper introduces the new method of Fourier synthesis pulse 
sequence design for systems showing dispersions in the 
parameters governing their dynamics. 

In this paper we focus on the Bloch equations with a linear 
one dimensional gradient, which arise in the context of NMR 
spectroscopy and NMR imaging applications. 










-G{t)s 


£u{t) 






My 




G(t)s 





~ev{t) 




My 


M, 




-eu{t) 


ev{t) 










(1) 



Here, M(s, e) = [M^ My M^Y is the state vector, u{t) e 
v{t) G 3?, and G{t) G 3? are controls and the parameters 
s G [0, 1] and e G [1 — (5, 1], S > are dispersion parameters 
which will be explained subsequently. Without loss of gener- 
ality, we will always normalize the initial state of the system 
([U to have unit norm, so that the system evolves on the unit 
sphere in three dimensions (Bloch sphere). A useful way to 
think about the Bloch equations is by imagining a two 
dimensional mesh of systems, each with a particular value 
of the pair (s, e). We are permitted to apply a single set of 
controls {u{t), v{t), G{t)) to the entire mesh of systems, and 
the controls should prepare the final state of each system as 
a desired function of the parameters (s, e) which govern the 
system dynamics. 

From a physics perspective, the system ([T]l corresponds 
to an ensemble of noninteracting spin-^ in a static mag- 
netic field Bq along the z axis and a transverse rf field 
{A{t) cos{ip{t)), A{t) sm{ip{t))) in the x-y plane. The state 
vector [Mx My MzY represents the coordinate of the unit 
vector in the direction of the net magnetization vector for 
the ensemble [4]. The controls u{t) and v{t) correspond to 
available rf fields we may apply to the ensemble of spins. 
The dispersion in the magnitude of the rf field applied to 
the sample is modeled by including a dispersion parameter 
e such that A{t) = eAo{t) with e G [1-5,1], 6 > 0. 
Thus, the maximum amplitude for the rf field (e — 1) 
corresponds to the maximum amplitude seen by any spin 
in the ensemble. Similarly, we consider a linear gradient 
G{t)s, where G{t) may be thought of as a control, and s 
represents the normalized spatial position of the spin system 
in the sample of interest. In ([T]), we work in units with the 
gyromagnetic ratio of the spins 7 = 1 . In this paper we give 
new design methods which scale polynomially that can be 
used to design pulse sequences for ([T]i which prepare the 
final state of the system as a function of the parameters s 
and e. 
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Fig. 1. A schematic depiction of tlie pulse sequence element in @. The 
pulse sequence element consists of six individual rotations which can be 
produced using the controls u{t) and v(t) as explained in the text. 



II. Design Method for re Inhomogeneity 

Considering only the Bloch equations with rf inhomogene- 
ity and no linear gradient {G{t) = 0), we can rewrite ([T]| in 
terms of the generators of rotation in three dimensions as 
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We will come back to the full version of the Bloch equations 
([T]i with both a linear gradient and rf inhomogeneity later in 
the paper. The problem is to design u{t) G and v{t) G 3? to 
effect some desired evolution. We now show how to construct 
controls to give a rotation of angle (/)(e) around the x axis 
or the y axis of the Bloch sphere. From these constructions, 
an arbitrary rotation on the Bloch sphere can be constructed 
using an Euler angle decomposition. 

A. Rotation About y axis 

In a time interval dt, we can use the controls to generate 
rotations exp{euQdtny) and exp^evodtflx) where uq and 
vq are constants to be specified. Using this idea, consider 
generating the rotation 



Uk = UikU2k 



(4) 



with 



Uik 

U2k 



exTp{-7rke^x) exp{-e(3k^y) exp{-KkeQx) (5) 
exp{Trkeil,x) exp{—ePk^y) exp(— Trfceilj;) (6) 



using the controls u and v. Using the relation 

exp(ai7a;) exp(/3f2y) exp(— afia;) = 

exp{P{cos{a)ny + sin{a)Qz)) (7) 

the matrices Uik and U2k may be rewritten as 

Uik — exp(ie/?fc(ilj^ cos(7rfce) — sin(7rfce))) (8) 
U2k = exp(-e/?fc(ilj^ cos(7rfce) + sin(7rfce))) (9) 



For small jSk, we can make the approximation 

UikU2k ~ exp{ePk coa{nke)fly) 



(10) 



In ( fTOb . we have expanded the exponentials in Uik and U2k 
to first order, performed the multiplication called for in (01, 
and then rewritten the product as ( fTOl l keeping terms to first 
order. In the case when /3k is too large for ( fTol i to represent 
a good approximation, we should choose a threshold value 
Po such that ( fTOl i represents a good approximation and so 
that 

l3k=n(3o (11) 
with n an integer. Defining 

Uio = exp(ie/3o(f^a cos(7rfce) — sin(7rfce))) (12) 

C^20 = exp(ie/3o(ilj^ cos(7rfce) + sin(7rfce))) (13) 

we can apply the total propagator 

[ C/ioC/20 ]" (14) 
« [exp(e/3o cos(7rfce)Sly)]"' (15) 
— exp{e(3k cos{Trke)Qy) (16) 

where we used the approximation (fTol i in dTsl ). More will be 
said about this approximation below. 

If we then think about making the incremental rotation Uk 
for many different values of k, we will get a net rotation 



[/ = cxp(e/3fc cos{Trke)i}y) 



(17) 



so long as we keep Pk sufficiently small to justify the 
approximation (fTOl i. The total propagator U for the Bloch 
equations can then be rewritten as 



U = cxp(e ^ Pk cos{'Kke)i}y) 



If we now choose the coefficients Pk so that 

Pk cos(7rfce) « 

k 



(18) 



(19) 



then we will have constructed a pulse sequence to approxi- 
mate a desired e dependent rotation around the y axis. Since 
e is bounded away from the origin, 0(e)/e is everywhere 
finite, and we can approximate it with a Fourier Series. 



1 ) Remark About the Approximation: Here we consider 
the error introduced by the approximation ( fTsT l. Define the 
error E{Z^ V) when a unitary matrix V is implemented 
instead of a desired unitary matrix Z by 



E{Z,V) =maxl|(Z-F)a;|| 

X 

With these identifications, in ([T4l > we have 

V = U10U20 

= I + — cos(nke)ny + MUn) 



(20) 



(21) 



and 



exp(/3o cos(7rfce)rij,) 

/ + — cos(7rfce)f7„ + Mzfn) 



(22) 



where Mi (n) and M2 (n) are matrices with finite entries and 
of maximum order . Notice this imphes that the difference 
{Z — V) is of order As defined previously, Pk — n(3o- 
A well known result (see for example [19]) says that the 
maximum error introduced by implementing (multiplying) 
the product of n of the V matrices instead of implementing 
n of the Z matrices is the sum of the individual errors. Thus, 
the total error -Etotai satisfies 



total 



< nE{Z, V) 



1 

n 



(23) 
(24) 



and thus, by making n sufficiently large, we may decrease the 
total error introduced in the above method to an arbitrarily 
small value. For the simulations done in this paper, we find 
a value /3o < 30 deg produces good results. 

B. Rotation About x axis 

An analogous derivation can be made for rotations about 
the X axis of the Bloch sphere. Replacing (|5]l and (|6]l with 

Uik = exp{~nkei}y) exp{^e(3k^x) exp{TTkeV,y)(25) 
U2k = exp(7rfceriy)exp(ie/3fcri^)exp(-7rfcefij^)(26) 

and following an analogous procedure leads to an approxi- 
mate net propagator 



U = cxp(e ^ /3fc cos(7rfce)ila;) 



(27) 



The coefficients f3k may be chosen to approximate 0(e)/e, 
and thus we can approximately produce a desired e depen- 
dent rotation about the x axis of the Bloch sphere. Since an 
arbitrary rotation on the Bloch sphere may be decomposed in 
terms of Euler angles, the methods presented can be used to 
approximately synthesize any evolution on the Bloch sphere. 



C. Choosing the Coefficients (3k 

Suppose we wish to design a pulse with a uniform net 
rotation angle of </> around either the x or y axis of the 
Bloch sphere using the previously discussed algorithm. We 
focus on the case of a uniform rotation (independent of e) 
because this is the most useful pulse sequence in NMR. It is 
straightforward to incorporate an e dependent rotation (f){e) 
into everything that follows. We face the problem of choosing 
Pk and k so that 

Pk cos(7rfce) « /(e), l~5<e<l (28) 



where 



/(e) = l-<5<e<l (29) 



Since we only have cos(7rfce) terms in the series, we first 
will extend /(e) to have even symmetry about e — 0. To do 
this, we define g(e) to be 

r /(e), l-5<e<l 
5(e) = < /(l--?), -(l-<5)<e<l-,5 (30) 
[ /(-e), -l<e<-(l-<^) 

and now consider choosing f3k and k so that 



k 



/3fcCos(7rfce) w g(e), 



-1< e < 1 



(31) 



A natural choice is to choose k as nonnegative integers, in 
which case (ik may be computed using the orthogonality 
relation 

f-i 

cos(7rfce) cos(7rfc'e)(ie = Jfefe' , 7^ (32) 
where 5kk' is the Kronecker delta. We find for the coefficients 

Pk — J cos{TTke)g{e) de, k ^ 



g{e) de, k = 



The number of terms kept in the series is decided by the pulse 
designer. A sufficient number of terms should be retained so 
that the error across the relevant range of e does not exceed 
some acceptable value. Figure |2] depicts an example design 
using a series with five terms. The region of interest for /(e) 
is 0.1 < e < 1. In this region, we see relatively small errors. 
We now give two examples to demonstrate the usefulness of 
the algorithm. 

D. Simulations 

1) pulse around y axis: Suppose we wish to design 
a ^ pulse around the y axis and we want to consider rf 
inhomogeneity in the range 0.1 < e < 1. Then we should 
consider 

/(e) = 0.1<e<l (33) 

Figure |3] shows the results of the designed pulse sequence 
acting on the initial state M(0) = [0 1]"^ while keeping 
five terms in the series expansion. We see that the resulting 
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Fig. 2. An example design to approximate over tlie range — 1 < e < 1 
for the rotation angle = t/Z. Five terms in the series expansion were 
retained. The relevant range of rf inhomogeneity is 0.1 < e < 1. 



1.2r 



Fig. 4. Results of applying u{t) = ^ for one unit of time to the initial state 
M(0) = [0 1]^. The system coiresponding to the value e = 1 for the 
dispersion parameter experiences a rotation about the y axis of the Bloch 
sphere, but systems corresponding to other values of e show deteriorated 
performance. 
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Fig. 3. Results of pulse sequence designed to produce uniform 7r/2 rotation 
about y axis. The sequence was applied to the initial state A/(0) = [0 1]^, 
and the plot shows the final state as a function of e after propagating the 
Bloch equations. 



Fig. 5. Results of pulse sequence designed to produce uniform tt rotation 
about X axis. The sequence was applied to the initial state A/(0) = [0 1 O]'^, 
and the plot shows the final state as a function of e after propagating the 
Bloch equations. 



pulse sequence reliably produces a net evolution exp(^rij,) 
across the entire range of e values. Figure |4] shows the results 
of applying u{t) = ^ for one unit of time to the system 
(|2]l. This approach corresponds to assuming the dispersion 
parameter e is fixed at a nominal value e = 1, so that 
every system sees the same control signals u{t) and v{t). 
Systems corresponding to values e 1 exhibit deteriorated 
performance as demonstrated in Figure ID 

2) TT pulse: As a second example, suppose we wish to 
design a tt pulse around the x axis and we want to consider 
rf inhomogeneity in the range 0.5 < e < 1. Then we should 
consider 

/(e) = 0.5<e<l (34) 



Figure |5] shows the results of the designed pulse sequence 
acting on the initial state M(0) = [0 1 0]^ while keeping 
nine terms in the series expansion. We see that the resulting 
pulse sequence reliably produces a net evolution exp(7rf7y) 
across the entire range of e values. 

It should be noted that although we consider design 
examples where we wish to produce a uniform rotation that is 
independent of the parameter e, the method presented in the 
paper can also be used to design control laws which prepare 
the final state as a function of the parameter e. We consider 
examples to produce a uniform rotation, independent of e, 
because this is the most useful application in NMR. 



III. Design Method for Position Dependent 
Rotations 

Now consider the Bloch equations with no rf inhomogene- 
ity and with a Hnear gradient 
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(35) 

where u{t) e 3?, v{t) e 5R, and G{t) e 3? are time 
dependent control ampHtudes we may specify, and s G [0, 1] 
can be thought of as a dispersion parameter As previously 
discussed, s represents the spatial position of the spin system 
in the sample under study. The goal is to engineer a control 
law that will effect a net position dependent rotation, so that 
the final state is prepared as a function of s. A common 
example in NMR imaging is a so-called slice selective 
sequence, whereby the controls should selectively perform 
a ^ rotation on some range of s values, while performing 
no net rotation on s values falling outside of that range. 
Using the controls, consider generating the evolution 



Uk = UikU: 



2k 



(36) 



with 



Uik = exp(7rfcsr2z)exp(-/3feriy)exp(-7rfcsr22)(37) 
U2k = exp{-TrksV,z)exp{-l3ki^y)exp{'Kksi}z){3S) 
The matrices Uik and U2k may be rewritten as 

Uik = exp(i/3fc(cos(7rfcs)f7,y - sin(7rfcs)ria;)) (39) 

U2k = ex-pi -f3k{cos{TTks) fly + sin(TTks)ilx)) (40) 

Again performing a first order analysis on the exponentials 
as was done in the previous section, we can make the 
approximation 



UikU2k ~ exp(/3fc cos(7rfcs)r2j^) 



(41) 



If we then think about making the rotation Uk for different 
values of k and f3k, we will get the net propagator 



[/ = JJ^ exp(/3fc cos(7rfcs)f7j^) 



(42) 



within the approximation previously discussed. The propa- 
gator ( l42t may be rewritten as 



U = exp(^ f3k cos{Trks)i}y) 



(43) 



Choosing k as the nonnegative integers, and choosing the (3k 
so that 

f3k cos(7rA:s) « 0(s) (44) 

k 

where (j){s) is the desired position s dependent rotation angle 
results in a net rotation around the y axis of the Bloch sphere 
with the desired dependence on the parameter s. 
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Fig. 6. Results of pulse sequence designed to produce uniform ^ rotation 
about y axis over the range 0.5 < s < 0.75. The sequence was applied to 
the initial state Af(0) = [001]"^, and the plot shows the final state as a 
function of s after propagating the Bloch equations. 



An analogous procedure may be used to generate an s 
dependent rotation around the x axis of the Bloch sphere. 
Replacing ^7} and (l38]l with 

Uik = exp{-TTksftz)exp{^(3k^x)exp{TTksilz)i4-5) 

U2k = exp{nksV,z)exp{^(3k^x)ew{-T^ksilz)(46) 

and following a similar procedure, we can approximately 
produce the total propagator 



U = exp(^ Pk cos{TTks)^j^ 

k 



(47) 



Choosing k as the nonnegative integers, and choosing the 
(3k appropriately results in a net rotation around the x axis 
of the Bloch sphere with the desired dependence on the 
parameter s. Since any rotation can be decomposed in terms 
of Euler angles, we may use the methods just discussed to 
approximately produce any position s dependent rotation on 
the Bloch sphere. 

A. Design Example 

As an example design using the procedure just discussed, 
consider a slice selective pulse sequence, where we wish to 
excite a certain range of s values while leaving systems with 
s values falling outside of that range unaffected at the end 
of the sequence. 



f , 0.5 < s < 0.75 

0, otherwise 



(48) 



Figure |6] shows the results of a pulse sequence designed using 
the procedure described in the text while keeping 30 terms 
in the series. The ripples appearing in Figure |6] result from 
the ripples in the approximation of the sharp slice selective 
profile using a Fourier Series. One method used to overcome 



this in practice is to allow for a ramp between the and ^ 
level on the slice. 

IV. Control Laws Involving Position and rf 
Inhomogeneity 

We now come back to the problem of considering the 
fuU version of the Bloch equations ([T]i including the two 
dispersion parameters s and e. Rewriting ([T]i in terms of the 
generators of rotation we have 



A choice of ki, k2, and (3k so that 



/3fc cos(7rfcis) cos(7rA;2e) ~ 



(55) 
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(49) 

where the state vector is now a function of both parameters 
s and e. The control task is to choose u{t) G 5R, v{t) G 
and G{t) G K to effect a desired rotation </)(s,e). Proceed- 
ing along the lines of the previous two sections, consider 
generating the propagators 

Uik — exp(7rfcisilz) exp(-e/3fer2j,) exp(— TrfcisSl^) 

= exp(-e/?fc(cos(7rfcis)rij, — sm{nkis)Qx)) (50) 



and 



U2k = exTp{—Trkisflz) exp{—e(3k^y) exp{TTkisQz) 

= exp{-ePk{cos{T:kis)fly + sin{nkis)Qx)) (51) 



Within a first order approximation for the exponentials, we 
have the approximate total propagator 

C^ifeC^2fc « exp(ie/3fcCos(7rfcis)rij^) (52) 

Building on this, we can produce the propagator 

Usk = exp(7rfc2ef2:c)[/ifeL/2fc exp(-7rfc2erix) 

w exp(-e/3fe cos(7rfcis)(cos(7rfc2e)riy — sin(7r/c2e)f^z)) 

Similarly, we can produce 

U4,k = e-Kp{-Trk2eflx)UikU2kexp{Trk2e^x) 

Ri cxp(-e/3fe cos(7rfcis)(cos(7rfc2e)f7j, + sm{Trk2e)flz)) 

so that we can approximately produce the total propagator 

Uk = Uj,kUik 

« exp(e/3fc cos(7rfcis) cos(7rfc2e)riy) (53) 

within the approximation for the exponentials. We can use 
the method previously discussed in the case when Pk is 
too large for the approximation to be valid. Producing the 
propagator Uk for different values of fci, fc2, and j3k results 
in the net propagator 

U= exp(e/?fe cos(7rA:is) cos(7rfc2e)rij;) (54) 

{kiM} 



where (/)(s, e) is the desired position s and rf inhomogeneity 
parameter e dependent rotation angle, results in an approxi- 
mate desired evolution for the Bloch equations ( |49] l. 

Analogous arguments show we may approximately pro- 
duce a rotation around the x axis of the Bloch sphere 

Uk — exTp{ef3k cos(7rfcis) cos(7rfc2e)f^2:) (56) 
and may thus approximately generate a net propagator 

U— exp(e/3fc cos(7rfcis) cos(7rfc2e)rix) (57) 

{kiM} 

and thus approximately produce a desired position s and 
rf inhomogeneity parameter e rotation around the x axis of 
the Bloch sphere. Since any rotation on the unit sphere can 
be decomposed in terms of Euler angles, an arbitrary (s, e) 
dependent rotation can be approximately produced using 
these methods. 

V. Conclusions 

In this paper we have provided new methods to design 
control laws for the Bloch equations when certain dispersion 
parameters are present in the system dynamics. These meth- 
ods are of utmost practical importance in the fields of NMR 
spectroscopy and NMR imaging, and can be implemented 
in many well known experiments immediately. The methods 
presented in the paper allow the design of a compensating 
control law (pulse sequence) that will compensate for disper- 
sions in the system dynamics while providing a clear tradeoff 
for the control law designer between total time required for 
the sequence and amplitude of the available controls. 
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